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TITLE OF THE INVENTION 

System and Method for the Automatic Extraction of Linear Features from Digital 
Imagery 

FIELD OF THE INVENTION 

[0001] The present invention relates to digital imagery. More 

specifically, the present invention is concerned with a system and a method for 
the automatic extraction of linear features from digital imagery. 

BACKGROUND OF THE INVENTION 

[0002] With the recent, rapid development of computer tools, digital 

image processing has become an efficient, economic, and convenient means 
for obtaining qualitative and quantitative information in different fields such as 
remote sensing, cartography, robotics, and materials. 

[0003] Indeed, digital image processing makes it possible to 

qualitatively describe images from various sources in terms of pattern 
recognition to identify and isolate contained objects. One of the major 
subclasses of contained objects present in digital images consists in linear 
features. Automatic detection of linear features from digital images plays an 
important role in pattern recognition and in digital image processing. 

[0004] A variety of techniques and many algorithms have emerged 

to automatically extract linear features from digital images. These techniques 
can be classified into two main categories: local methods which are based on 
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local operators such as mobile kernel, and global methods which focus on 
mathematical transformations, such as Hough transforms (see [1]). 

[0005] Local methods for automatically extracting linear features 

basically exploit local variations of pixel intensity in a small neighborhood by 
calculating the gradients in small, limited-size windows in the image, e.g., 3x3 
pixels or 5 x 5 pixels (see [2] [3] [4] [5] [6] [7] [8] [9]). 

[0006] A number of researchers have examined mathematical 

morphology as a means of extracting linear features (see [10][11][12][13][14]). 
The problems pertaining to this technique arise from the number of human 
decisions required to reconnect and rebuild line segments, which increases 
processing time. Multi-dimensional line detection is the other technique for 
detecting for linear features that collects different spectral information for the 
same scene and may highlight different parts of lines (see [15]). The first stage 
of this method for obtaining the combined images requires several 
transformations of multiple original bands. Human intervention is needed to 
select the best-combined image. Another approach to linear-feature extraction 
involves knowledge-based systems, which need more information than a 
simple digital image for line extraction (see [16][17][18]). 

[0007] These local methods generally remain inefficient because 

they fail to have a global view of the linear features in a digital image. One 
problem common to all of these methods is that the resulting extracted line 
images contain a fair amount of noise, while the detected lines are often 
incomplete and geometrically shifted. These difficulties are magnified by 
intersections between different linear features and linear features that display 
some curvature (see [19][20][1]). In addition, these methods turn in 



exceptionally long processing time when extracting features from large images 
(see [16]). 

[0008] Of the local methods, LINDA (Linear-featured Network 

Detection and Analysis) system (see [2][7][27]), based on the profile intensity 
analysis of the pixel line, is the most recent and the most representative 
method. With all of basic disadvantages of the local methods, the LINDA 
system is far from being operational with respect to systematically and 
automatically processing a large set of images. 

[0009] The Radon transform and its derivative, the Hough transform, 

are the most frequently used approaches as global methods for detecting linear 
features (see [21][22][23][24]). In principle, a straight line from the input image 
is transformed into a digital peak (a light or dark pixel, compared to its 
neighborhood) in the transformed plane. In this case, it is easier to detect a 
peak in the transformed plane than a straight line in the input image. There are 
three basic limitations for these methods that unfortunately restrict their 
applications and their utility in practice. 

[0010] First, the Hough and Radon transform-based methods are 

limited solely to extracting straight lines (see [21]). Therefore, linear features 
that span the entire image but display some curvature may not produce 
suitable peaks or troughs in the transform plane. This restriction is linked 
directly to the basic definition of the method. 

[001 1] Secondly, there are false peaks or virtual maxima (peaks that 

do not represent the real line on the input image) in the transformed plane. 
These false peaks considerably decrease the quality of the results by 
increasing the error of commission. When the line density in the input image is 
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high, eliminating the false peaks from transformed plane seriously limits these 
methods. 

[0012] The last methodological inconvenience of these approaches 

is the ignorance of the nature of the detected lines. Indeed, since the intensity 
integration in the transformation process is performed over the entire length of 
the image, the process can have difficulty detecting line segments that are 
significantly shorter than the image dimensions (see [21][25][26]). Neither can it 
provide information about the positions of the endpoints of these shorter line 
segments or line length. 

[0013] It is to be noted that the present description refers to a 

number of documents indicated in brackets, that are listed hereinbelow, the 
content of these documents is herein incorporated by reference in their entirety. 

OBJECTS OF THE INVENTION 

[0014] An object of the present invention is therefore to provide an 

improved system and method for the automatic extraction of linear features 
from digital imagery. 

SUMMARY OF THE INVENTION 

[0015] More specifically, in accordance with the present invention, 

there is provided a method for the extraction of linear features from digital 
imagery comprising the steps of: 

providing a digital image; 

providing a multi-layer database 

initializing a parameter domain; 
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successively applying the Radon transform on each position in 
the parameter domain; 

for each position in the parameter domain: 

- finding the coordinates of the nearest pixel in the digital image; 

- determining the numerical value of the found pixel; and 

- populating the layers of the provided database in accordance 
with the determined numerical value of the found pixel; 

analysing the data of the database; and 

generating an output image by restoring the lines detected in the 
analysing step. 

[0016] According to another aspect of the present invention, there is 

provided a system for the extraction of linear features from a digital image and 
for generating a corresponding output image comprising: 

a controller provided with an input designed to receive a digital 
image and an output designed to provide a corresponding output image; said 
controller containing a multi-layer database; said controller being so configured 
as to: 

initialize the discrete parameter domain; 

successively apply a Radon transform on each position in the 
parameter domain; 

for each position in the parameter domain: 

- find the coordinates of the nearest pixel in the digital image; 

- determine the numerical value of the found pixel; and 

- populate the layers of the provided database in accordance 
with the determined numerical value of the found pixel; 

analyse the data of the database; and 

generate an output image by restoring the lines detected in the 
analysing step. 



[0017] It is to be noted that the expression "linear features" used 

herein and in the appended claims is to be construed as referring both to 
curved lines as well as straight lines. 

[0018] Other objects, advantages and features of the present 

invention will become more apparent upon reading of the following non- 
restrictive description of preferred embodiments thereof, given by way of 
example only with reference to the accompanying drawings. 

BRIEF DESCRIPTION OF THE DRAWINGS 
[0019] In the appended drawings: 

[0020] Figures 1A and 1B illustrate a flowchart of an embodiment of 

the method of the present invention; 

[0021] Figure 2 is an image simulated to showcase the method's 

ability to extract all lines without consideration for the slope line; extract the 
lines having a shorter length of the image dimension; and locate extracted lines 
accurately; 

[0022] Figure 3 is an image obtained by applying the implemented 

algorithm to the image of Figure 2; 

[0023] Figure 4 is an image simulated in order to demonstrate the 

method's performance in extracting curved as well as straight lines; 

[0024] Figure 5 is an image produced by applying the method of 

Figures 1A and 1B to the image of Figure 4; 
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[0025] Figure 6 is a SPOT satellite image of a part of a road map; 

[0026] Figure 7 is an image produced by applying the method of 

Figures 1 A and 1B to the image of Figure 6; 

[0027] Figure 8 is a scanned image of a part of a fingerprint; 

[0028] Figure 9 is an image produced by applying the method of 

Figures 1A and 1B to the image of Figure 8; and 

[0029] Figure 10 schematically illustrates a system according to an 



|H embodiment of the present invention. 

U DESCRIPTION OF THE PREFERRED EMBODIMENT 

Mi 

G [0030] The proposed methodology generally derives from the Radon 

j[7 transform, yet has none of its basic disadvantages mentioned hereinabove. 
This methodology also builds on a solid mathematical base that demonstrates 
the conceptual feasibility of the new method mathematically. An algorithm has 
also been developed and implemented in order to achieve a concrete 
realization of the innovative method in a real context. 

[0031] Generally stated, the proposed methodology is based on the 

tridimensional Radon transformation and on the generation of a multi-layer 
database having separate layers containing different types of information to 
allow the reconstruction of an output image where the linear features of the 
input image are restored. 
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[0032] In following subsections, the mathematical base as well as 

general procedures for an algorithmic development of an embodiment of the 
method and system of the present invention will be presented in details. 

Mathematical Development 

Continuous function 

[0033] Letg(x,y) be a continuous signal of the continuous variables x 

and y and let t, denote a r|-dimensional parameter vector defined as: 

(1) 

[0034] where ^ spans the parameter domain. 

[0035] For a two-dimensional continuous function, the parameter 

vector can be defined as: 

$ = (^,§3) 1 = 3 (2) 

or: 

^ = (a,P,y) (3) 

hence: 

g(x,y) = <K*,y;£) = y-ax 2 -px - y (4) 



[0036] 



The Radon transform g(§) of function g(x 9 y) is defined as: 
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§0s) = f f g( x > y) 5 (<f>(x, y;^)dxdy (5) 
[0037] where 8(.) denotes the Dirac delta function. 

[0038] Using the definition (5) of Radon transform, curves expressed 

by the next parameter form can be detected: 

4>(x,y;5) = 0 (6) 

[0039] Substituting the relation (4) into equation (5), we find: 



g(a,p,y)= ££g(x,y)8(y-ax 2 ~px-y)dxdy (7) 
[0040] where y = <|>(x;£) represent the transformation curve. 

[0041] Using proprieties of delta function, equation (7) becomes: 

g(ct,P,y) = ££g(x,y)5(y-ax 2 -Px-y)dxdy- £g(x 5 ax 2 +Px + y) dx (8) 

[0042] At this stage, the transformation curve is defined as a 

polynomial function of second degree using the delta function as below: 

g(x,y) = 8(y-<Kx£*)) (9) 
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[0043] 



where the parameter vector of this curve is: 



r=(a*,p\y) 



(10) 



[0044] 



Using equations (4), (9) and (10), we find: 



g(x,y) = 5(y-(oc*x 2 + p*x + y*)) 



(11) 



[0045] 



According to definition of the Radon transform: 



g(Q= f8(<Kx;S)-<Kx;5*))dx 

J- 00 

=f.i 5(x ' Xi) 



SMsS at**;?*) 



= 1 



dx dx 



1 



dx 



dx dx 



(12) 



[0046] Substitution of parameter vectors (3) and (10) into the 

transform definition (12) gives: 



g(a,p,y) = £$((ax 2 + px + y) - (a * x 2 + p * x + y*)) dx 



8(x-x ; ) 



g(ax 2 +px + y) g(a*x 2 +p*x + y*) 
dx dx 



dx (13) 
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1 



d(ocx 2 + px + 


d(a*x 2 +p*x+y*) 


dx 


dx 



[0047] Finally, we find: 

g(a,j3,y)=eo for a = a* et j3 = /3* et y = y* 

(14) 

0 -< g(a, f3,y)<co far other cases 

[0048] Equation (14) shows that the Radon transform of a 

polynomial function of second degree will give an infinite value (a peak) in the 
parameter domain, when parameters of this function are precisely determined. 
In addition, this equation confirms that the parameter domain will contain some 
other non-zero values because of parameters values which have some 
similarity with real parameter values of the polynomial function in question. 

[0049] In the light of the mathematical demonstration presented 

above, tree important issues are concluded: 

i) The Radon transform of a two-dimensional plane containing a 
polynomial function of second degree is a Euclidean tree-dimensional 
space; 

ii) this polynomial function is mapped into a unique point in the tree- 
dimensional space; 

iii) position of this point in the tree-dimensional space determines the 
polynomial function parameters. 
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.pass. 



[0050] The following subsection provides a discrete approximation to 

the Radon transform for a polynomial function of second degree. This discrete 
approximation is the key of the algorithmic development of the polynomial 
function in question. 

Discrete Function 

[0051] Discrete Radon transform of a discrete function (|>(m;l,k,h) is 

defined as: 



g(l,k,h) = £g(m,<Km;U,h))= £g(m,n(m;l,k,h)) (15) 

m=0 m=0 

where : 

<|>(m;l,k,h) = n = l"»( x - ^^MX^-y^)' 

L A y 

and 



(16) 



f X m =X m in+ mAx m = 0,l,--,M 

[y n = y nm + nAy n =0,1,--,N 



(17) 



[0052] where M and N are the discrete plane dimensions and Ax, Ay 

are sampling intervals of this plane. 

[0053] For a polynomial function of second degree expressed as 

y - ex 2 - px - x = 0 , the Radon transform is: 

M-l , 

g(l, k, h) = Ax ^ g(x m , s,x^ + p k x m + x h ) (18) 
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with: 



6i = + lAs 
Pk = Pmm + kAp 



1 =0,1,-,L 
k = 0,l,—,K 
h = 0,l,--,H 



(19) 



[0054] where L, K and H represent dimensions of Radon domain 

and As, Ap and Ai are sampling intervals of this domain. 



[0055] 

becomes: 



Using relation (19), the function n(m;l,k,h) of relation (16) 



n(m;l,k,h): 



g l X m+Pk X m+Th-yn 

Ay 



(20) 



[0056] 

we find: 



With substituting the relation (17) in the previous equation, 



n(m;l,k,h) = 



Si( x min + ™Ax) 2 + p^x^ + mAx) + x h - y n 
Ay 



= n - 



(21) 



[0057] 



For simplifying, the previous equation can be expressed as: 



n* = am 2 + pm + y 



(22) 
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where: 



a = 



SjAx 



P Ay 

2 

y = 8 l X min + Pk X min + T h ~ 

Ay 



(23) 



[0058] According to the equation (19), once the Radon domain 

parameters has been initialized, using equations (21), (22) and (23) this domain 
is determined for all of two-dimensional plane points (X,Y). 

Algorithmic development 

[0059] The concrete achievement of the proposed method is 

schematized in a flowchart presented in Figures 1A and 1B. This flowchart 
illustrates also the instructions for an algorithmic development of the proposed 
technology. 

[0060] As will easily be understood by one skilled in the art, the 

algorithm presented hereinbelow aims at building a multi-layer database. A 
first layer of the database includes the coordinates of pixels having a zero value 
and a second layer includes the coordinates of the pixels having a non-zero 
value. A third layer of the database includes the values of the pixels having a 
non-zero value while the fourth layer of the database includes the value of the 
pixel having a zero value. Finally, a fifth layer of the database includes the sum 
of the values of the non-zero value pixels. 
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[0061] In a first step 10, the digital image 12 and the minimun and 

maximum lengths (in pixels) of lines to be detected 14 are supplied to design 
the discrete parameter domain. The parameter vector of the transformation is 
defined in this step as described in equations (1 ), (2) and (3) hereinabove. 

^fe^ 2 ,4 3 -^,-^) (1) 
5 = tfU 2 ,S 3 ) n = 3 (2) 

S = (a,P,Y) (3) 



[0062] A second step 16 consists in the initialization of the numerical 

values of discrete parameters s mi n, As, p min , Ap, i min , Ai, L, K et H, as described 
in equation (19) and of x m , y n , Ax, Ay, as described in equation (17) (M and N 
corresponding to the dimension of digital input image). 



I y n = Yn™ + nAy 



m = 0,l,---,M 
n =0,1,—,N 



(17) 



°1 ~ °min 
Pk =Pmi 



+ lAs 
+ kAp 



1 =0,1,-,L 
k = 0,l,---,K 
h = 0,l,— ,H 



(19) 



[0063] In a third step 18, a first loop is started. This first loop 

intends, for every position in the parameter domain, to determine every 
possible combination of si, p k and t h as described in equation (19) by 
incrementing tree values of I, k and h (0 < I < L; 0 < k< K; and 0 < h < H). This 
first loop includes a second loop as will be discussed hereinbelow. 
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[0064] For each position in the parameter domain, a computation 

step 20 is performed. More specifically, the computation step 20 consists in a 
computation of line/curve deterministic parameters. In this step, the 
deterministic parameters a, (3 and y are computed by using equation (23). 



a = 



P = 



SjAx 2 

( 2X mm £ l+Pk) Ax 

Ay 

s i x L + Pk x mfa+' c h-y n 

Ay 



(23) 



[0065] In step 22, a second loop is started. This second loop 

generally consists in the determination of corresponding pixel's coordinates and 
in the update of the various layers of the database described hereinabove. In 
this second loop, for each value of m (0<m<M), a new parameter n* is 
computed by using equation (21). 



n(m; l,k,h) = 



SlfonA, + 1"A X ) 2 + Pk( X m ,n + mA X ) + T h ~ y n 

Ay 



(21) 



[0066] In the next step 24, the determined coordinates (m, n*) of the 

previous step are used to find the nearest neighbour approximation pixel P(m, 
n) within the input image. 
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[0067] Then, step 26 determines if the value of the pixel is zero. 

More specifically, this is a verification step for determining the numerical value 
of the pixel P(m, n) found in step 24. 

[0068] If the value of the pixel P(m, n) is zero, a number of zero pixel 

counter is incremented (step 28), the pixel's coordinates are saved in the first 
layer of the database (step 30) and the value of the pixel is saved in the fourth 
layer of the database (step 32). The counter incremented in step 28 is 
generally used to determine the gap between pixels that have a non-zero 
value. 

[0069] It is to be noted that when the input image used includes 

shades of grey, step 26 may be viewed as a comparison between the 
numerical value of the pixel and a predetermined threshold to determine if this 
numerical value may be considered as being zero. It is in this case that it is 
advantageous to save the value of the zero-value pixel in the fourth layer of the 
database. 

[0070] The system is then returned to the start of the second loop 

(step 22) for the next position in the parameter domain. 

[0071] If the value of the pixel P(m, n) is not zero, a number of non- 

zero pixel counter is incremented (step 34), the pixel's coordinates m and n are 
saved in the second layer of the database (step 36) and the numerical value of 
the pixel is saved in the third layer of the database (step 38). The counter 
incremented in step 34 is generally used to determine the length of the lines 
defined by the pixels. 
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[0072] The numerical pixel value is also added to the previous pixel 

values found (step 40), and the accumulated pixel value is stored in the fifth 
layer of the database (step 42). The fifth layer of the database is generally used 
to determine the length of the line defined by the pixels. 

[0073] When the second and first loops are completed, i.e. when the 

previously discussed steps have been performed for every position in the 
parameter domain, the method calls for a third loop that is started in step 44. 

[0074] Generally stated, the function of the third loop is to recall all 

the numerical values of non-zero pixels already stored in the third layer of the 
database in step 38. The third loop is therefore done for each pixel stored in 
the third layer of the database. 

[0075] In step 46, the numerical value of the pixel is compared to 

predetermined value relating to the length of the lines to be detected inputted in 
step 10. More specifically, the zero pixel and non-zero pixel counters are used 
to determined the length of the line represented by the third layer cell is within 
the range determined by the minimum and maximum length of lines to be 
detected 14. 

[0076] If it is not the case, the method calls for a return to the third 

loop if it has not reached its end. 

[0077] If it is the case, step 48 recalls the stored coordinates of the 

corresponding pixels in the third layer of the database. Then, in the input 
image, the numerical values of these pixels will be verified to find endpoints of 
the corresponding line. 
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[0078] Then, in step 50, by using the stored deterministic 

parameters of the line and the line endpoints that are found in the precedent 
step, the detected line will be restored in the output image. 

[0079] Finally, all the information about extracted lines, such as 

endpoints, lengths, etc. are stored (step 52) in a new database which is used to 
generate an output image including the linear features of the corresponding 
input image. 

Features and Operation of the Invention 

[0080] A global method has therefore been developed so as to be 

capable of: 

- detecting and extracting the lines from a digital image with any curvature; 

- discarding the virtual maxima on the Hough or on the Radon transformed 
plane; 

- locating every line segment in the input image (endpoints and length). 
Demonstration and Discussion 

[0081] The image of Figure 2 was simulated to showcase the 

method's ability to: 

i) extract all lines without consideration for the slope line; (for example line 
100) 

ii) extract the lines having a shorter length than that of the image 
dimension (for example line 102); and 

iii) locate extracted lines accurately. 
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[0082] Figure 3 is an image obtained by applying the above- 

described method to the image of Figure 2. The absolute similarity of Figure 3 
with Figure 2 demonstrates clearly the accuracy of three abilities of the method 
as described above. 

[0083] The image of Figure 4 was simulated in order to demonstrate 

the method's performance in extracting curved (for example line 104 as well as 
straight lines (for example line 106). 

[0084] The image illustrated in Figure 5 was produced by applying 

the above-described algorithm to the image of Figure 4. Again, the absolute 
similarity of Figure 5 with Figure 4 demonstrates clearly the accuracy of the 
method as described above. 

[0085] As a visual demonstration, a urban road network (see Figure 

6) obtained from a SPOT satellite (from CNES (Centre National d'Etudes 
Spatiales) has been selected as an application. 

[0086] The image of Figure 6 is a part of SPOT satellite image of a 

road map. The image of Figure 7 was produced by applying the above- 
described algorithm to the image of Figure 6. 

[0087] A comparison of the images of Figures 6 and 7 reveals that 

the restrained roads in the input image have been extracted in detail and with 
high geometric precision. 
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[0088] Turning now briefly to Figures 8 and 9 of the appended 

drawings, another application of the method and system of the present 
invention will be described. Figure 8 illustrates a portion of a scanned 
fingerprint. On the other hand, Figure 9 illustrates the output image when the 
input image is Figure 8. As will be apparent to one skilled in the art, the output 
image of Figure 9 could be used in a fingerprint identification system. 

[0089] Finally, Figure 10 of the appended drawings schematically 

illustrates a system 200 for the extraction of linear features from an input 
image. The system 200 includes a controller 202, for example in the form of a 
personal computer (not shown) that is so designed to run a computer program 
embodying the method described hereinabove. The controller 202 therefore 
includes the multi-layer database. 

[0090] The system 200 also includes an input device 204, for 

example in the form of a scanner (not shown) or of a link to another computer 
containing the input image (not shown) accessible through a computer network. 

[0091] An output device 206 is also provided to output the image. 

This output device may be in the form of a printer (not shown), a computer 
monitor (not shown) or of a link to another computer configured to receive the 
output image (not shown) accessible through a computer network. 

[0092] The method and system of the present invention are 

advantageous since they provide for the automatically extraction of line 
features and provide accurate, acceptable results. Furthermore, the method of 
the present invention is advantageous since it can extract curved lines of any 
length with a high degree of accuracy. Our method identifies every line or line 
segment, labeling its length, its curvature, and endpoint coordinates in a 
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database. Moreover, the database can be easily integrated into any geographic 
or spatial information system. 

Examples of Fields of Applications 

[0093] The method and system of the present invention is of broad 

interest and make a general contribution to the field of pattern recognition and 
digital image processing, regardless of the nature or source of the digital 
image. It could be used in a vast number of fields and application, of which the 
following are a sampling. 

[0094] In the field of remotely sensed imagery, the method and 

system of the present invention could be used for the extraction and update of 
road maps; for planimetric cartography, geology, and hydrology; for mapping 
power-lines, pipelines, and railways; and for tracking cycle trails, ORV trails, 
and the like, for example. 

[0095] In the field of medical imagery, the method and system of the 

present invention could be used for the detection of arteries, blood vessels, 
tumors, for example. 

[0096] In the field of imagery from other sources, the method and 

system of the present invention could be used for electrical circuits and the like. 

[0097] As will be apparent to one skilled in the art, even though the 

method has been presented herein with a single multi-layer database, it could 
also be implemented with a number of databases each representing a layer 



23 

described hereinabove. The expression "database layer 5 ' is therefore to be 
construed as also meaning a separate database. 

[0098] Furthermore, the method and system of the present invention 

may be used on binary and "grayscale" images. 

[0099] Although the present invention has been described 

hereinabove by way of preferred embodiments thereof, it can be modified, 
without departing from the spirit and nature of the subject invention as defined 
in the appended claims. 
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